function out = simProjection(model,shocks)
% Simulating the following variables (in deviation from steady state)

numSim       = size(shocks,2);
ny           = model.ny;
nx           = model.nx;
nx1          = model.nx1;
mx           = model.mx;    %number of pure endogenous states
myx          = model.myx;   %number of controls that are states;
nx1          = model.nx1;   %number of endogenous states

ySim         = zeros(ny,numSim);
xSim         = zeros(nx,numSim);
heteroShocks = reshape(model.heteroShocks,model.ne,1);

if model.pruningOn == 1
    xf_cu    = zeros(nx,1);
    xs_cu    = zeros(nx,1);
    xrd_cu   = zeros(nx,1);
    xAllSim  = zeros(nx+(model.appOrder-1)*nx1,numSim);
    for s = 1:numSim
        if model.appOrder == 1
            ySim(:,s) = model.g0 + model.gx*xf_cu;
            xSim(:,s) = xf_cu;
            xAllSim(1:nx,s) = xf_cu;
        elseif model.appOrder == 2
            AA_cu     = kron3(xf_cu,xf_cu);
            ySim(:,s) = model.g0 + model.gx*(xf_cu+xs_cu) + model.Gxx*AA_cu;
            xSim(:,s) = xf_cu+xs_cu;
            xAllSim(1:nx+nx1,s) = [xf_cu;xs_cu(1:nx1,1)];
        elseif model.appOrder == 3
            AA_cu     = kron3(xf_cu,xf_cu);
            BB_cu     = kron3(xf_cu,xs_cu);
            ySim(:,s) = model.g0 + model.gx*(xf_cu+xs_cu+xrd_cu) + model.Gxx*(AA_cu + 2*BB_cu) ...
                        + model.Gxxx*kron3(xf_cu,AA_cu);
            xSim(:,s) = xf_cu+xs_cu+xrd_cu;
            xAllSim(1:nx+2*nx1,s) = [xf_cu;xs_cu(1:nx1,1);xrd_cu(1:nx1,1)];
        end
        
        % Updating the states (dev from steay state)
        xf_cup = model.h0 + model.hx*xf_cu;
        
        xf_cup(nx1+1:nx,1) = xf_cup(nx1+1:nx,1) +...
                (2./(exp(-heteroShocks(:).*xf_cu(nx1+1:nx,1))+1)).*model.eta(nx1+1:nx,:)*shocks(:,s);
                         
        xf_cu = xf_cup;
        if model.appOrder > 1
            xs_cup = model.hx*xs_cu + model.Hxx*AA_cu;
            xs_cu  = xs_cup;
        end
        if model.appOrder > 2
            xrd_cup = model.hx*xrd_cu + model.Hxx*(2*BB_cu) ...
                + model.Hxxx*kron3(xf_cu,AA_cu);
            xrd_cu  = xrd_cup;
        end 
    end
else
    x_cu = zeros(nx,1);
    for s = 1:numSim
        xSim(:,s)   = x_cu;
        if model.appOrder == 1
            ySim(:,s) = model.g0 + model.gx*x_cu;
        elseif model.appOrder == 2
            ySim(:,s) = model.g0 + model.gx*x_cu + model.Gxx*kron3(x_cu,x_cu);
        elseif model.appOrder == 3
            AA_cu = kron3(x_cu,x_cu);
            ySim(:,s) = model.g0 + model.gx*x_cu + model.Gxx*AA_cu + model.Gxxx*kron3(AA_cu,x_cu);
        end
        % Updating the states
        if model.mx > 0
            % update endogenous states (not lagged controls) - in dev from steady state 
            x_cup(1:model.mx,1) = updateEndoStates(model,x_cu);
        end
        x_cup(mx+1:mx+myx,1) = ySim(1:myx,s) - model.gSS(1:myx,1); %- in dev from steady state
        x_cup(nx1+1:nx,1)    = model.h0(nx1+1:nx,1) + model.hx(nx1+1:nx,:)*x_cu ...
                             + (2./(exp(-heteroShocks(:).*x_cu(nx1+1:nx,1))+1)).*model.eta(nx1+1:nx,:)*shocks(:,s);
        
        x_cu  = x_cup;
    end
    
end

%% Saving output
out.y_cu  = ySim;
out.x_cu  = xSim;
if model.pruningOn == 1
    out.xAll_cu = xAllSim;
end

end